


%
clear; clear all;
data = readmatrix('Orignal_time_delay_45deg_100Hz.txt'); 
%% Read Data
% Assume your text file data format is: first column time (seconds), second column time delay (as)
%data = load('your_data_file.txt'); % Replace 'your_data_file.txt' with your filename
time = data(:, 1);    % First column: time (s)
delay = data(:, 2);    % Second column: time delay (as)

%% Data Preprocessing
% Remove mean to analyze fluctuations around the mean. Skip if analyzing absolute fluctuations.
delay_detrended = delay - mean(delay);

delay_to_use = delay_detrended;
Fs = 1 / (time(2) - time(1)); % Calculate sampling frequency (Hz), assuming uniform sampling interval
%% Method 1: Using pwelch function (Recommended, uses averaging and windowing for smoother, more stable results)
% Set parameters
window_length = min(256, floor(length(delay_to_use)/2)); % Window length, typically a fraction of data length
overlap_ratio = 0.5; % Overlap ratio, e.g., 50%
nfft = max(1024, 2^nextpow2(length(delay_to_use))); % FFT points, typically a power of 2

%% Method: Using periodogram function (Direct method)
[psd_per, freq_per] = periodogram(delay_to_use, [], nfft, Fs);
psd_comparison = [freq_per, psd_per, 10*log10(psd_per)];
header = 'Frequency(Hz)    PSD_Welch(linear)    PSD_Welch(dB)';
save('PSD_45deg_100Hz.txt', 'psd_comparison', '-ascii', '-tabs');



%
clear; clear all;
data = readmatrix('Orignal_time_delay_45deg_0.1Hz.txt'); 
%% Read Data
% Assume your text file data format is: first column time (seconds), second column time delay (as)
%data = load('your_data_file.txt'); % Replace 'your_data_file.txt' with your filename
time = data(:, 1);    % First column: time (s)
delay = data(:, 2);    % Second column: time delay (as)

%% Data Preprocessing
% Remove mean to analyze fluctuations around the mean. Skip if analyzing absolute fluctuations.
delay_detrended = delay - mean(delay);

delay_to_use = delay_detrended;
Fs = 1 / (time(2) - time(1)); % Calculate sampling frequency (Hz), assuming uniform sampling interval
%% Method 1: Using pwelch function (Recommended, uses averaging and windowing for smoother, more stable results)
% Set parameters
window_length = min(256, floor(length(delay_to_use)/2)); % Window length, typically a fraction of data length
overlap_ratio = 0.5; % Overlap ratio, e.g., 50%
nfft = max(1024, 2^nextpow2(length(delay_to_use))); % FFT points, typically a power of 2

%% Method: Using periodogram function (Direct method)
[psd_per, freq_per] = periodogram(delay_to_use, [], nfft, Fs);
psd_comparison = [freq_per, psd_per, 10*log10(psd_per)];
header = 'Frequency(Hz)    PSD_Welch(linear)    PSD_Welch(dB)';
save('PSD_45deg_0.1Hz.txt', 'psd_comparison', '-ascii', '-tabs');



%
clear; clear all;
data = readmatrix('Orignal_time_delay_01deg_100Hz.txt'); 
%% Read Data
% Assume your text file data format is: first column time (seconds), second column time delay (as)
%data = load('your_data_file.txt'); % Replace 'your_data_file.txt' with your filename
time = data(:, 1);    % First column: time (s)
delay = data(:, 2);    % Second column: time delay (as)

%% Data Preprocessing
% Remove mean to analyze fluctuations around the mean. Skip if analyzing absolute fluctuations.
delay_detrended = delay - mean(delay);

delay_to_use = delay_detrended;
Fs = 1 / (time(2) - time(1)); % Calculate sampling frequency (Hz), assuming uniform sampling interval
%% Method 1: Using pwelch function (Recommended, uses averaging and windowing for smoother, more stable results)
% Set parameters
window_length = min(256, floor(length(delay_to_use)/2)); % Window length, typically a fraction of data length
overlap_ratio = 0.5; % Overlap ratio, e.g., 50%
nfft = max(1024, 2^nextpow2(length(delay_to_use))); % FFT points, typically a power of 2

%% Method: Using periodogram function (Direct method)
[psd_per, freq_per] = periodogram(delay_to_use, [], nfft, Fs);
psd_comparison = [freq_per, psd_per, 10*log10(psd_per)];
header = 'Frequency(Hz)    PSD_Welch(linear)    PSD_Welch(dB)';
save('PSD_01deg_100Hz.txt', 'psd_comparison', '-ascii', '-tabs');



%
clear; clear all;
data = readmatrix('Orignal_time_delay_01deg_0.1Hz.txt'); 
%% Read Data
% Assume your text file data format is: first column time (seconds), second column time delay (as)
%data = load('your_data_file.txt'); % Replace 'your_data_file.txt' with your filename
time = data(:, 1);    % First column: time (s)
delay = data(:, 2);    % Second column: time delay (as)

%% Data Preprocessing
% Remove mean to analyze fluctuations around the mean. Skip if analyzing absolute fluctuations.
delay_detrended = delay - mean(delay);

delay_to_use = delay_detrended;
Fs = 1 / (time(2) - time(1)); % Calculate sampling frequency (Hz), assuming uniform sampling interval
%% Method 1: Using pwelch function (Recommended, uses averaging and windowing for smoother, more stable results)
% Set parameters
window_length = min(256, floor(length(delay_to_use)/2)); % Window length, typically a fraction of data length
overlap_ratio = 0.5; % Overlap ratio, e.g., 50%
nfft = max(1024, 2^nextpow2(length(delay_to_use))); % FFT points, typically a power of 2

%% Method: Using periodogram function (Direct method)
[psd_per, freq_per] = periodogram(delay_to_use, [], nfft, Fs);
psd_comparison = [freq_per, psd_per, 10*log10(psd_per)];
header = 'Frequency(Hz)    PSD_Welch(linear)    PSD_Welch(dB)';
save('PSD_01deg_0.1Hz.txt', 'psd_comparison', '-ascii', '-tabs');





%% Visualize Results
figure;
subplot(2,1,1);
plot(time, delay);
xlabel('Time (s)');
ylabel('Time Delay (as)');
title('Original Time Delay Signal');
grid on;

subplot(2,1,2);
%plot(freq_welch, 10*log10(psd_welch), 'b-', 'LineWidth', 1.5); hold on;
plot(log10(freq_per), 10*log10(psd_per), 'r--', 'LineWidth', 1);
%legend('Welch Method', 'Periodogram Method');
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)'); % Note: dB is decibel, a logarithmic scale
title('Power Spectral Density (PSD) of Time Delay');
grid on;
